clc,clear;
close all;
% 种植面积
data=xlsread('2023种植面积data.xlsx');
data(isnan(data))=0;
% 预测销售量
sale=readmatrix('预测销售量.xlsx');
sale=sale(:,3);
% 亩产量
chag=readmatrix('亩产量数据表格.xlsx');
chag=chag(:,3:end);
% 销售单价
r=readmatrix('销售单价数据.xlsx');
r=r(:,3:end);
% 种植成本
p=readmatrix('种植成本数据.xlsx');
p=p(:,3:end);
% 地区最大面积
S=xlsread('附件1.xlsx'); 
S=S([1:54,27:54],1);
% 最低种植面积
ms=readmatrix('最低种植面积.xlsx');
ms=ms(:,3:end);
%参数
y=7;%年份
m=82;%地区类型
n=41;%作物类型
maxi=20000;%最大迭代次数
T0=1000;%初始温度
a=0.999;%降温率

x0=data;%初始数据
f0=profit(x0,chag,sale,p,r);
for k=1:y
    if k==1
        x=x0; %更新初始数据
    end
    T=T0;
    count=1;
    %进行优化模型求解(模拟退火)
    for i=1:maxi
        %产出新解
        newx=x;
        %每次对于一类地区中的作物进行一定的扰动
        %(1-26)地区可(1-15)作物,(27-54)可(17-34),产物16只能27-34
        %(55-62)可35-37,(63-78)可38-41,(79-82)可(17-34)
        %分为6类
        op=randi([1,6]);%随机的类
        if op==1
            i0=[1,26];j0=[1,15];
        elseif op==2
            i0=[27,54];j0=[17,34];
        elseif op==3
            i0=[27,34];j0=16;
        elseif op==4
            i0=[55,62];j0=[35,37];
        elseif op==5
            i0=[63,78];j0=[38,41];
        else
            i0=[79,82];j0=[17,34];
        end
       pl=0.5;
      if pl>rand
            if j0==16
                io=randi(i0);
                newx(io,j0)=rand*newx(io,j0)+ms(io,j0);
            else
                io=randi(i0);
                jo=randi(j0,3,1);
                for lk=1:3
                newx(io,jo(lk))=rand*newx(io,jo(lk))+ms(io,jo(lk));
                end
            end
      else
             %交换作物
             if op~=3
                 for hi=min(j0):max(j0)
                     for ki=min(i0):max(i0)
                          if newx(ki,hi)>0
                              msk=newx(ki,hi);
                          else
                              continue;
                          end
                          jo=randi(j0);
                          while jo==hi
                              jo=randi(j0);
                          end
                          mst=newx(ki,jo);
                          newx(ki,jo)=msk;
                          newx(ki,hi)=mst;
                     end
                 end
             end
       end
        %对于豆类提供更容易的概率选中
        po=0.5;
        if po>rand
            jj=1:5;
            ji=[1,5];
            for ki=1:26
                if sum(newx(ki,jj))==0
                    s=randi(ji);
                    newx(ki,s)=ms(ki,s);
                    for kj=6:n
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
            ji=[17,19];
            jj=17:19;
            for ki=27:50
                if sum(newx(ki,jj))==0
                    s=randi(ji);
                    newx(ki,s)=ms(ki,s);
                    for kj=[1:16,20:n]
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
            for ki=51:54
                if sum(newx(ki,jj))==0||sum(newx(ki+28,jj))
                    s=randi(ji);
                    if rand>0.5
                        newx(ki,s)=ms(ki,s);
                    else
                        newx(ki+28,s)=ms(ki+28,s);
                    end
                     for kj=[1:16,20:n]
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
        end
        %检查是否满足约束条件
        if k==1
        [newx,flag]=yuesu(newx,S,ms,k,x0);
        elseif k==2
           [newx,flag]=yuesu(newx,S,ms,k,x0,y1);
        elseif k==3
            [newx,flag]=yuesu(newx,S,ms,k,y1,y2);
        elseif k==4
            [newx,flag]=yuesu(newx,S,ms,k,y2,y3);
        elseif k==5
            [newx,flag]=yuesu(newx,S,ms,k,y3,y4);
        elseif k==6
            [newx,flag]=yuesu(newx,S,ms,k,y4,y5);
        elseif k==7
            [newx,flag]=yuesu(newx,S,ms,k,y5,y6);  
        end
        if flag==0
         continue;
        elseif flag==1&&count==1
            %计算适应值并更新最优值
            newf=profit(newx,chag,sale,p,r);
            %设定初始条件
            bestx=newx;
            bestf=newf;
            f=newf;
            x=newx;
            count=0;
        end
        %计算适应值并更新最优值
        newf=profit(newx,chag,sale,p,r);
        delta=newf-f;
        if delta>0||exp(-delta/T)<rand
            x=newx;
            f=newf;
            % 更新全局最优
             if newf>bestf
                 bestx=newx;
                 bestf=newf;
             end
        end
        T=T*a;
    end
    disp(bestf)
    if k==1
        y1=bestx;
    elseif k==2
        y2=bestx;
    elseif k==3
        y3=bestx;
    elseif k==4
        y4=bestx;
    elseif k==5
        y5=bestx;
    elseif k==6
        y6=bestx;
    elseif k==7
        y7=bestx;
    end
    x=x0;
    bestx=x;
    bestf=f0;
    f=f0;
end
function f=profit(x,ch,y,p,r)
   f=0;
   for i=1:41
        if sum(x(:,i).*ch(:,i))-y(i)>0
            t=1;
        else
            t=0;
        end
        f=f+sum(x(:,i).*ch(:,i).*r(:,i)-x(:,i).*p(:,i))-sum(t*r(:,i).*(sum(x(:,i).*ch(:,i))-y(i)));
   end
end
function [x,flag]=yuesu(x,S,ms,k,yx1,yx2)
   [m,n]=size(x);
   flag=1;
   for i=1:m
           %种植面积约束
           if x(i,:)>S(i)
               flag=0;
               return
           end
       for j=1:n
           %正向约束
           if x(i,j)<0
               flag=0;
               return
           end
       end
   end 
   %最小种植面积与地区约束
   for i=1:26  %第一大类
      for j=1:n
          if j>15
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif j<=15&&x(i,j)>0
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
           if  j<=15
              if k==1
                  if yx1(i,j)*x(i,j)>0
                    x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end  
              end
          end
      end
   end
   for i=27:34  %第二大类
      for j=1:n
          if j<16||j>34
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&& j>=16&&j<=34
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
          if  j>=16&&j<=34
              if k==1
                  if yx1(i,j)*x(i,j)>0
                    x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end  
              end
          end
      end
   end
   for i=35:54  %第三大类
      for j=1:n
          if j<17&&j>34
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&& j>=17&&j<=34
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
           if  j>=17&&j<=34
              if k==1
                  if yx1(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end  
              end
          end
      end
   end
   for i=55:62  %第四大类
      for j=1:n
          if j<35&&j>37
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&& j>=35&&j<=37
               if x(i,j)<ms(i,j)
               flag=0;
                return
               end
          end
          if  j>=35&&j<=37
              if k==1
                  if yx1(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end  
              end
          end
      end
   end
   for i=63:78  %第五大类
      for j=1:n
          if j<38&&j>41
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&& j>=38&&j<=41
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
           if j>=38&&j<=41
              if k==1
                  if yx1(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                         x(i,j)=0.9*x(i,j);
                  end  
              end
          end
      end
   end
   for i=79:82  %第六大类
      for j=1:n
          if j<17&&j>34
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&& j>=17&&j<=34
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
          if  j>=17&&j<=34
              if k==1
                  if yx1(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                        x(i,j)=0.9*x(i,j);
                  end  
              end
          end
      end
   end
   %三年豆类约束
   j=1:5;
   for i=1:26
       if k>1
           if sum(yx1(i,j))+sum(yx2(i,j))+sum(x(i,j))==0
               flag=0;
               return
           end
       end
   end
   j=17:19;
   for i=27:50
       if k>1
           if sum(yx1(i,j))+sum(yx2(i,j))+sum(x(i,j))==0
               flag=0;
               return
           end
       end
   end
   for i=51:54
       if k>1
           if sum(yx1(i,j))+sum(yx2(i,j))+sum(x(i,j))==0&&sum(yx1(i+28,j))+sum(yx2(i+28,j))+sum(x(i+28,j))==0
               flag=0;
               return
           end
       end
   end
end


